Section: New Results

Computational Statistical Physics

Participants : Matthew Dobson, Claude Le Bris, Frédéric Legoll, Tony Lelièvre, Francis Nier, Grigorios Pavliotis, Mathias Rousset, Gabriel Stoltz.

The extremely broad field of molecular dynamics is a domain where the MICMAC project-team, originally more involved in the quantum chemistry side, has invested a lot of efforts in the recent years. Molecular dynamics may also be termed computational statistical physics since the main aim is to numerically estimate average properties of materials as given by the laws of statistical physics. The project-team studies both deterministic and probabilistic techniques used in the field. One of the main difficulty is related to the metastable features of the generated trajectories: the system remains trapped over very long times in metastable states, which means that very long trajectories need to be generated in order to obtain macroscopically relevant quantities. This is related to the fact that the timescale at the microscopic level is much smaller than the timescale at the macroscopic level. In [66] , we propose a summary of the mathematical approaches to quantify metastability, and which appear to be useful to analyze the numerical methods used in molecular dynamics.

Free Energy calculations

For large molecular systems, the information of the whole configuration space may be summarized in a few coordinates of interest, called reaction coordinates. An important problem in chemistry or biology is to compute the effective energy felt by those reaction coordinates, called free energy.

In the article [42] , Tony Lelièvre, Mathias Rousset and Gabriel Stoltz study the application of constrained Langevin dynamics to the computation of free energy differences, by thermodynamic integration techniques and fluctuation relation (à la Jarzynski).

One interest of free energy computation techniques is that they appear to be useful in other fields, like in computational statistics where multimodal measures are also frequently encountered, so that standard Markov Chain Monte Carlo appraoches also suffer from metastability.

For example, in  [25] , Nicolas Chopin (CREST, ENSAE), T. Lelièvre and G. Stoltz explore the application of the Adaptive Biasing Force method to Bayesian inference. This sampling method belongs to the general class of adaptive importance sampling strategies which use the free energy along a chosen reaction coordinate as a bias. Such algorithms are very helpful to enhance the sampling properties of Markov Chain Monte Carlo algorithms, when the dynamic is metastable.

In [58] , G. Fort (Telecom Paris), B. Jourdain (CERMICS), E. Kuhn (INRA), T. Lelièvre and G. Stoltz have considered the Wang-Landau algorithm. The authorshave proved that the Wang-Landau algorithm converges with an associated central limit theorem, and have provided an analysis of the efficiency of the algorithm in a metastable situation.

Convergence to equilibrium

An important question for the analysis of sampling techniques is the rate of convergence to equilibrium for stochastic trajectories.

In [65] , F. Nier, T. Lelièvre and G. Pavliotis study the interest of using non-reversible stochastic dynamics to enhance the rate of convergence to equilibrium, compared to reversible dynamics. A well posed optimization problem is obtained and solved in the case of a linear drift for the overdamped Langevin dynamics.

Metropolis Hastings algorithms

A classical sampling tool used in molecular dynamics and in computational statistics is the Metropolis-Hastings algorithm. There has been a lot of work (see G. Roberts et al.) to study how the variance of the proposal should scale with the dimension of the problem, in order to optimize the sampling procedure. Most of these works assume that (i) the target probability is the product of $n$ one dimensional laws and that (ii) the Markov chain starts at equilibrium.

In the two works [60] , [59] , T. Lelièvre and his co-authors have generalized these results when the initial distribution is not the target probability. The diffusive limit in the latte case is solution to a stochastic differential equation nonlinear in the sense of McKean. They have discussed practical counterparts in order to optimize the variance of the proposal distribution to accelerate convergence to equilibrium. The analysis confirms the interest of the constant acceptance rate strategy (with acceptance rate between $1/4$ and $1/3$) first suggested in the works of G. Roberts et al., at least for the Random Walk Metropolis algorithm.

Thermodynamic limit

The quasicontinuum method is an approach to couple an atomistic model with a coarse-grained approximation in order to compute the states of a crystalline lattice at a reduced computational cost compared to a full atomistic simulation.

In that framework, the team has addressed questions related to the finite temperature modeling of atomistic systems and derivation of coarse-grained descriptions, such as canonical averages of observables depending only on a few variables. In the one-dimensional setting, an efficient strategy that bypasses the simulation of the whole system had been proposed in 2010. We refer to [47] for a recent review. In collaboration with X. Blanc (Université Pierre et Marie Curie), F. Legoll has extended this strategy to the so-called membrane setting in [16] .

When the temperature is small, a perturbation approach can be used to compute the canonical averages of these observables depending only on a few variables, at first order with respect to temperature. In collaboration with E. Tadmor, W. K. Kim, L. Dupuy and R. Miller, F. Legoll has analyzed such an approach in [46] . The numerical tests reported there show the efficiency of the approach, as long as the temperature is indeed small.

Sampling trajectories

There exist a lot of methods to sample efficiently Boltzmann-Gibbs distributions. The situation is much more intricated as far as the sampling of trajectories (and especially metastable trajectories) is concerned.

Following a numerical observation in a previous work on the sampling of reactive trajectories by a multilevel splitting algorithm, F. Cérou (Inria Rennes), A. Guyader (Inria Rennes), T. Lelièvre and F. Malrieu (Université de Rennes) study theoretically in [56] the distribution of the lengths of these trajectories, using large deviation techniques.

In [37] , C. Le Bris and T. Lelièvre together with M. Luskin and D. Perez from Los Alamos National Laboratoy provide a mathematical analysis of the parallel replica algorithm, which has been proposed by A. Voter in 1997 to simulate very efficiently metastable trajectories. This work opens a lot of perspectives, by using a generic tool (the quasi stationary distribution) to make a link between a continuous state space dynamics (Langevin dynamics) and a discrete state space dynamics (kinetic Monte Carlo models).

In a work in progress, T. Lelièvre and F. Nier have studied the quasi-stationnary distribution in relation for an overdamped Langevin process in a bounded domain. In the small temperature limit and by making the connection with boundary Witten Laplacians, they are able to compute accurately the spatial exit law along the boundary and non perturbative accurate formulas when the potential is changed inside the domain.

Effective dynamics

For a given molecular system, and a given reaction coordinate $\xi :{ℝ}^n\mapsto ℝ$, the free energy completely describes the statistics of $\xi \left(X\right)$ when $X\in {ℝ}^n$ is distributed according to the Gibbs measure. On the other hand, obtaining a correct description of the dynamics along $\xi$ is complicated.

F. Legoll and T. Lelièvre have introduced and analyzed some years ago a strategy to define a coarse-grained dynamics that approximates $\xi \left(X_t\right)$, when the state of the system $X_t$ evolves according to the overdamped Langevin equation (which is ergodic for the Gibbs measure). We refer to [47] for a recent review. The aim was to get a coarse-grained description giving access to some dynamical quantities (and not only equilibrium quantities). Together with G. Samaey (KU Leuven), they have recently studied how to use this coarse-grained description, accurate when the time scale separation is asymptotically large, to somewhat precondition the dynamics of the actual system in cases when the time scale separation is not large. For that purpose, they have used the parareal framework, to iteratively correct the sequential coarse-grained trajectory by fine scale trajectories performed in parallel. The main difficulty is that the two models (the reference one and the coarse-grained one) do not act on the same variable: the reference model evolves all the variables, whereas the coarse-grained model only evolves the slow variables. As shown in [63] in a simplified context (that of singularly perturbed ODEs), the precise coupling between both models should be done carefully.

The above study is concerned with models with continuous state spaces. S. Lahbabi and F. Legoll have studied in [61] a related question in the framework of kinetic Monte Carlo models, where the state space is discrete. For some models involving some slow and some fast variables, the effective dynamics of the slow component has been identified, and a complete proof of convergence proposed.

Hamiltonian dynamics

Constant energy averages are often computed as long time limits of time averages along a typical trajectory of the Hamiltonian dynamics. One difficulty of such a computation is the presence of several time scales in the dynamics: the frequencies of some motions are very high (e.g. for the atomistic bond vibrations), while those of other motions are much smaller. This problem has been addressed in a two-fold manner.

Fast phenomena are often only relevant through their mean effect on the slow phenomena, and their precise description is not needed. Consequently, there is a need for time integration algorithms that take into account these fast phenomena only in an averaged way, and for which the time step is not restricted by the highest frequencies. In [29] , M. Dobson, C. Le Bris, and F. Legoll have developed integrators for Hamiltonian systems with high frequencies. The integrators were derived using homogenization techniques applied to the Hamilton-Jacobi PDE associated to the Hamiltonian ODE. This work extends previous works of the team. The proposed algorithms can now handle the case when the (unique) fast frequency depends on the slow degrees of freedom, or when there are several fast constant frequencies.

Another track to simulate the system for longer times is to resort to parallel computations. An algorithm in that vein is the parareal in time algorithm. It is based on a decomposition of the time interval into subintervals, and on a predictor-corrector strategy, where the propagations over each subinterval for the corrector stage are concurrently performed on the processors. Using a symmetrization procedure and/or a (possibly also symmetric) projection step, C. Le Bris and F. Legoll, in collaboration with X. Dai and Y. Maday, have introduced several variants of the original plain parareal in time algorithm [28] . These variants, compatible with the geometric structure of the exact dynamics, are better adapted to the Hamiltonian context.

Nonequilibrium systems

The efficient simulation of molecular systems is known to be a much more complicated problem when the system is subjected to a non-conservative external forcing than when the system experiences conservative forces. Together with the sampling of metastable dynamics mentioned above, these are the two major research focus in molecular dynamics of the project-team.

Nonequilibrium molecular dynamics simulations can be used to compute the constitutive relation between the strain rate and stress tensor in complex fluids. This is fulfilled simulating molecular systems subject to a steady, non-zero macroscopic flow at a given temperature. Starting from a bath model, M. Dobson, F. Legoll, T. Lelièvre, and G. Stoltz have derived a Langevin-type dynamics for a heavy particle in a non-zero background flow [57] . The resulting dynamics, which is theoretically obtained when a unique large particle is considered, is numerically observed to also perform well when a system of many interacting particles within shear flow is considered.

Let us also mention that the article on the computation of the viscosity of fluids using steady state nonequilibrium dynamics with an external nongradient bulk forcing, in the framework of the PhD of Rémi Joubaud, has also been published [34] . In addition, the study by G. Stoltz and C. Bernardin on thermal transport in one-dimensional chains of oscillators whose kinetic and potential energy functions are the same, has been accepted and is now published [13] .